.. currentmodule:: brian

.. index::
   pair: example usage; IdentityConnection
   pair: example usage; NeuronGroup
   pair: example usage; run
   pair: example usage; exp
   pair: example usage; SpikeCounter
   pair: example usage; StateMonitor

.. _example-frompapers-computing with neural synchrony-duration selectivity_Fig1D_duration_selectivity:

Example: Fig1D_duration_selectivity (frompapers/computing with neural synchrony/duration selectivity)
=====================================================================================================

Brette R (2012). Computing with neural synchrony. PLoS Comp Biol. 8(6): e1002561. doi:10.1371/journal.pcbi.1002561
------------------------------------------------------------------------------------------------------------------
Figure 1C,D. Duration selectivity.
(takes about 1 min)

Caption (Fig. 1C,D).
A postsynaptic neuron receives inputs from A and B. It is
more likely to fire when the stimulus in the synchrony receptive field of
A and B.

::

    from brian import *
    
    # Parameters and equations of the rebound neurons
    Vt=-55*mV
    Vr=-70*mV
    El=-35*mV
    EK=-90*mV
    Va=Vr
    ka=5*mV
    gmax2=2
    tau=20*ms
    ginh_max=5.
    tauK2=100*ms
    N=10000 # number of neurons (= different durations)
    rest_time=1*second # initial time
    tmin=rest_time-20*ms #for plots
    tmax=rest_time+600*ms
    
    eqs='''
    dv/dt=(El-v+(gmax*gK+gmax2*gK2+ginh)*(EK-v))/tau : volt
    dgK/dt=(gKinf-gK)/tauK : 1 # IKLT
    dgK2/dt=-gK2/tauK2 : 1 # Delayed rectifier
    gKinf=1./(1+exp((Va-v)/ka)) : 1
    duration : second
    ginh = ginh_max*((t>rest_time) & (t<(rest_time+duration))) : 1
    tauK : ms
    gmax : 1
    theta : volt # threshold
    '''
    
    neurons=NeuronGroup(2*N,model=eqs,threshold='v>theta',reset='v=Vr;gK2=1')
    neurons.v=Vr
    neurons.theta=Vt
    neurons.gK=1./(1+exp((Va-El)/ka))
    # Neuron A, duplicated to simulate multiple input durations simultaneously
    neuronsA=neurons[:N]
    neuronsA.tauK=400*ms
    neuronsA.gmax=1
    neuronsA.theta=-55*mV
    neuronsA.duration=linspace(100*ms,1*second,N)
    # Neuron B, duplicated to simulate multiple input durations simultaneously
    neuronsB=neurons[N:]
    neuronsB.tauK=100*ms
    neuronsB.gmax=1.5
    neuronsB.theta=-54*mV
    neuronsB.duration=linspace(100*ms,1*second,N)
    
    # Noisy coincidence detectors
    tau_cd=5*ms
    tau_n=tau_cd
    sigma=0.2 # noise s.d. in units of the threshold
    eqs_post='''
    dv/dt=(n-v)/tau_cd : 1
    dn/dt=-n/tau_n+sigma*(2/tau_n)**.5*xi : 1
    '''
    postneurons=NeuronGroup(N,model=eqs_post,threshold=1,reset=0)
    CA=IdentityConnection(neuronsA,postneurons,'v',weight=0.5)
    CB=IdentityConnection(neuronsB,postneurons,'v',weight=0.5)
    
    spikes=SpikeCounter(postneurons)
    M=StateMonitor(postneurons,'v',record=N/3)
    
    run(rest_time+1.1*second,report='text')
    
    # Figure
    subplot(121) # Fig. 1C, example trace
    plot(M.times/ms,M[N/3],'k')
    xlim(1350,1500)
    ylim(-.3,1)
    xlabel('Time (ms)')
    ylabel('V')
    
    subplot(122) # Fig. 1D, duration tuning curve
    count=spikes.count
    # Smooth the tuning curve
    window=200
    rate=zeros(len(count)-window)
    for i in range(0,len(count)-window):
        rate[i]=mean(count[i:i+window])
    plot((neuronsA.duration[window/2:-window/2]/ms)[::10],rate[::10],'k')
    xlim(0,1000)
    ylim(0,0.5)
    xlabel('Duration (ms)')
    ylabel('Spiking probability')
    show()